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ANNUAL  RETORT 


ARTIFICIAL  STIMULATION  UF  EARTHQUAKES 
NEW  MEXICO  INSTITUTE  UF  MINING  AND  TECHNOLOGY 
Merle  Hanson 

Contract  No.  Fiji*  62Q-67-C-0113 

The  possibility  of  predictably  creating  artificial  earthquakes  through 
fluid  infection  into  geologic  structures  is  the  purpose  of  this  investigation, 
Seismic  activity  occurring  in  the  Denver,  Colorado,  region  which  may  or  may  not 
have  been  caused  by  injecting  fluids  into  the  Rocky  Mountain  Arsenal  waste 
injection  well  has  raised  the  question.  An  understanding  of  the  phenomena  was 
undertaken  through  analytic  means.  Obviously  this  work  has  a  direct  application 
for  the  objectives  of  underground  nuclear  test  control. 

This  report  contains  a  brief  summary  of  the  reasons  why  a  stable  two  dim¬ 
ensional  elastic  dynamics  numerical  calculations!,  technique  had  to  be  developed. 
In  addition,  the  technique  is  described  briefly  in  the  general  section  of  the 
report  and  in  more  detail  in  the  appendix.  The  results  discuss  some  simulations 
of  brittle  fractures  propagating  through  an  isotropic  elastic  media.  Emphasis 
is  placed  on  the  fact  that  che  sirrilatei  fractures,  using  sample  failure 
criteria,  propagate  at  speeds  greater  than  those  normally  predicted  in  the 
literature.  These  simulations  exhibit  the  fact  that  the  description  of  dynamic 
fracture  phenomenon  is  more  complex  than  those  used. 

Code  solution  comparison  with  analytic  solutions  in  seismic  problems 
indicate  that  the  method  is  sound.  Further  development  of  fracture  models  is 
required  to  understand  the  effect  of  high  pressure  fluids  in  creating  or  driving 
fractures  in  geologic  strata. 

A  photo  stress  meter  was  purchased  during  the  contract  year  to  use  for 
correlation  purposes  with  the  analytic  model. 


ABSTRACT 


'"'A  discussion  of  a  two  dimensional  elastic  dynamics  technique 
developed  to  simulate  seismic  disturbances  as  a  result  of  various 
Input  or  grid  created  disturbances.  The  method  exhibits  calculation*! 
stability  over  other  techniques  using  quadrilateral  grids.  A  quieter 
mesh  results  from  which  more  information  can  be  extracted  from  the 
calculation.  The  results  section  discusses  simulations  of  brittle 
fractures  propagating  through  the  isotropic  elastic  media.  Enphasis 
is  placed  on  the  fact  that  the  simulated  cracks,  using  simple  failure 
criteria,  propagate  at  speeds  greater  than  normally  predicted  in  the 
literature.  These  simulations  uxhibit  that  the  description  of  dynamic 
failure  is  more  complex  thin  those  used , 


INTRODUCTION  AND  SUMURY 


This  report  summarizes  the  work  performed  during  the  first  contract 
year  in  the  research  on  Artificial  Stimulation  of  Earthquakes,  Contract 
No.  F44620-67-C-0113,  at  the  New  Mexico  Institute  of  Mining  and  Technology. 
Initially,  the  work  was  directed  toward  applying  the  elastic  dynamics 
description  in  the  existing  two-dimensional  langrangian  tensor  hydrocodes 
to  the  seismic  problems  requiring  solution  for  this  analysis.  These 
methods,  when  an  appropriate  damping  was  used,  proved  to  be  caleulationally 
stable;  however,  the  danping  terra  was  not  adequate  for  the  low-power 
level  of  disturDances  calculated  for  the  seismic  problems;  it  tended  to 
spread  a  disturbance  over  large  areas  of  the  grid  in  a  nonphysical  manner. 
To  oyercoot>  this  problem,  equations  describing  a  Cosserat  continuum  were 
studied,  but  it  was  not  physically  realistic  to  apply  the  higher  order 
derivatives  in  terms  of  differences  to  the  IargrangiAn  merit.  Two 
alternatives  evident  at  this  time  for  applying  the  dynamic  finite  element 
technique  to  the  seismic  problem  were  (a)  to  use  computers  with  larger 
active  storage,  permitting  the  calculation  of  more  essh  intersections  and 
thereby  reducing  the  zone  size  so  that  the  spreading  of  the  disturbance, 
wnile  affecting  the  same  number  of  zones,  would  cover  a  much  smaller  part 
of  the  continuum  o-‘  (b)  to  try  to  develop  another  description  that  would 

include  additional  constraints  in  the  equations  describing  strain,  thereby 
superseding  the  requirement  for  the  damping  constraints  needed  in  sob© 
of  the  existing  techniques.  For  obvious  reasons,  the  decision  was  made 
to  develop  a  calculation*!  technique  suitable  for  the  simulation  of  seismic 
type  of  problems.  This  technique  was  developed  and  is  described  briefly 
in  th&»  appendix.  Correlation  with  problems  having  analytic  solutions  has 
not  shown  serious  differences  between  the  code  solutions  and  the  analytic 
solutions.  During  the  latter  part  of  the  contract  year,  this  code  has 
beer,  applied  to  problems  of  hydraulically  driven  cracks  and  shear  cracks, 
’lynamic  simulations  of  crack  phenomena  using  the  code  do  not  compa/e  in 
crack  velocity  with  generally  accepted  puolisned  results  indicating  that 
dynamic  failure  of  isotropic  elastic  media  is  quite  complicated.  A  lit¬ 
erature  search  is  in  progress  for  information  that  will  help  provide  an 
understanding  of  the  cracking  phenomena. 
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The  possibility  of  predictably  creating  artificial  earthquakes 
through  fluid  injection  into  geologic  structures,  is  the  purpose  of 
this  investigation.  Earthquake  activity  in  the  Denver,  Colorado,  area 
during  ard  if  ter  waste-fluid  injeotion  in  the  Rocky  Mountain  Arsenal 
injection  well  indicated  that  artificial  triggering  of  earthquake  activity 
injecting  high-pressuiu  fluid  into  geologic  strata  aught  oe  possible. 

We  do  not  intend  to  imply  that  the  injection  of  fluids  at  the  RM1  well 
caused  the  seismic  activity  in  the  Denver  area;  however,  it  is  certainly 
conceivable.  An  understanding  of  the  phenomena  was  undertaken  through 
analytic  means  by  simulating  stress  fields  dynamically  with  a  two-dimensional, 
finite  element  technique.  Although  the  analysis  has  not  been  completed, 
the  approach  seems  sound  and  in  addition  to  solutions  of  dynamic  release 
of  stress  fields  in  rock  needed  for  this  task,  the  technique  should  prove 
valuable  for  other  seismological  studies.  Obviously,  this  work  has  a 
direct  application  for  underground  nuclear  test  control,  since  a  signif¬ 
icant  background  of  seismic  noise  created  either  intentionally  or 
unintentionally  could  mask  underground  nuclear  shots.  In  addition,  small 
seismic  disturbances  could  be  misunderstood  for  underground  nuclear  shots 
or  vice  versa. 
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GENERAL 


A  s+.anaard  method  for  handling  two-dimensional  elastic -plastic 
formulations  in  numerical  form  has  been  published  by  Wilkins  [l]  . 

By  his  method,  the  strain  rates  within  a  Iangrangian  mesh  ar® 
calculated  from  the  velocities  and  co-ordinates  of  the  mesh  poiits, 
then  Integrated  to  obtain  the  strain.  This  method  does  not  provide 
the  necessary  constraints  for  a  complete  description  of  the  strain  ard, 

t 

hence,  the  stress  field.  For  this  reason,  the  calculations  are 
absolutely  unstable,  so  the  method  is  sometimes  stabilized  by  use  of 
what  is  called  a  rotational  damping  term  [2].  The  rotational  damping 
constraints  are  obtained  by  calculating  pseudo  pressures  from  the  rates 
of  rotation  of  line  segments  on  both  sides  of  the  mass  point.  The  pseudo 
pressures,  which  are  proportional  to  difference  of  the  rates  of  rotation 
of  opposing  line  segments  of  a  mesh  point,  are  then  applied  to  oppose 
the  motion  of  the  mesh  point  to  reduce  the  difference  in  the  upposlng 
rates  of  rotation.  Hiis  method  appears  to  work  well  where  the  event  to 
be  simulated  occurs  in  short  time  periods  compared  to  the  rotation  or 
distortion  of  the  mesh  caused  by  the  lack  of  constraint  in  the  strain 
description. 

This  technique  for  constructing  difference  equations  in  iAngrangian 
form,  when  applied  to  seismic  problems,  spreads  a  disturbance  In  a  non¬ 
physical  manner  over  large  areas  of  the  grid,  making  it  extremely  complex. 

The  rotational  damping  term,  which  calculates  pseudc  pressures  of  the  order 
of  me'jnitude  of  the  calculated  stress  level  In  the  grid  near  the  disturb¬ 
ance,  caused  the  spreading.  Dimensions  of  the  extended  source  were  large 

I 

y  compared  to  a  physically  realistic  detector  dimension.  Without  the 

rotational  damping  terms,  the  grid  distorted  nonphysically  to  the  point 
where  little  information  could  be  extracted  from  the  calculations.  This 
left  two  alternatives;  either  to  attempt  to  build  another  formulation 
where  the  constraints  required  for  stability  were  contained  in  the  numerical 
algorithms  or  to  attempt  to  use  the  existing  formulations  on  larger  computers, 
which  would  allow  the  zoning  to  be  sufficiently  fine  that  the  rotational 
damping  terms,  while  affecting  the  same  number  of  zones,  would  cover  only 
small  parts  of  the  total  mesh  and  would  not  seriously  affect  the  results 
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of  the  calculation.  Reduction  of  tne  zone  size  in  the  continuum 
description  would  also  reduce  the  calcuiational  time  step  per  cycle, 
thereby  increasing  the  number  of  cycles  to  be  calculated  for  a  given 
problem.  Increasing  the  number  of  cycles  along  with  thi  greater  number 
of  zones  required  to  describe  the  continuum  would  make  the  computer 
running  time  prohibitively  long.  In  addition,  one  wants  to  describe  the 
largest  continuum  possible  with  the  conputer  equipment  available,  sine*; 
a  nonphysical  boundary  will  always  create  reflections  from  incident  waves 
in  an  elastic  madium.  Since  it  is  not  usually  possible  to  i  unulate  the 
complete  continuum  to  the  detail  required  with  existing  computer  equip¬ 
ment,  only  small  sections  can  be  considered,  thereby  introducing  non¬ 
physical  boundaries  where  the  section  of  the  continuum  terminates,  foe 
might  absorb  the  energy  at  the  boundary  from  an  incident  wave  with  a 
given  wave  number,  but  in  a  two-dimensional  description  of  elastic  media, 
waves  of  different  numbers  nearly  always  exist.  For  example,  consider 
the  irrotational  and  equivolume  waves  in  a  Hookean  material. 

Because  we  wanted  to  simulate  the  largest  segment  possible  of  the 
continuum  to  obtain  the  most  information  before  reflected  waves  from  the 
boundaries  affect  the  media  of  interest,  we  directed  work  toward  develop¬ 
ment  of  a  formulation  that  could  adequately  perform  the  calculation. 

This  resulted  in  a  new  two-dimensional  elastic  dynamics  code,  briefly 
described  in  the  appendix. 

Correlation  of  this  technique  with  problems  having  analytic  solutions 
has  not  shown  serious  differences  between  the  code  solutions  and  the  an¬ 
alytic  solutions.  Th©  following  observations  have  been  made  on  a 
computer  simulation  of  Lamb's  Problem  (a  source -detect  or  problem  on  an 
elastic  isotropic  half  space): 

1.  P,  S,  and  Rayleigh  waves  exist  and  were  calculated  at  nearly 
the  expected  velocities. 

2.  A  point  source  seems  to  be  extended  only  over  the>  zone  size 
of  the  Langrangian  mesh. 

3.  The  calculations  were  stable  without  iny  damping,  either  tensor 
or  rotational. 
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4.  The  direction  of  motion  resulting  from  the  P,  S,  And  Rayleigh 
waves  agrees  with  the  expected  motions. 

Other  observations  have  been  made  during  checKs  and  correlation 
work  with  the  coda.  The  method  calculates  the  correct  stress  per  zone 
in  both  the  shear  and  normal  cases,  dynamically,  this  code  simulates 
slightly  reduced,  but  predicted,  wave  speeds.  According  to  a  stability 
analysis,  all  finite  elemental  calculational  techniques  wilx  show  a 
reduced  wave  speed.  The  amount  a  wave  is  slowed  is  a  function  of  the 
zone  size,  input  disturbance  rise  time,  and  the  calculational  time  step. 
Dispersion  errors  can  lie  expected  if  an  input  disturbance  has  a  rise 
time  shorter  than  the  fastest  wave  traversal  time  for  a  Langrangian  zone. 
Other  checks  are  described  In  reference^ 

The  computational  code  can  handle  up  to  2400  mesh  points  on  an 
IBM  360  Model  44  computer.  Computational  speed  with  the  code  in  its 
present  form  is  about  4500  zone-cycles  a  second  on  this  computer. 

The  two-dimensional  technique  is  now  being  applied  to  the  phenomenon 
of  hydraulically  driven  fractures  In  brittle  elastic  isotropic  media.  In 
addition,  strata  under  high  shear  stress  with  preferred  lines  of  failure 
are  being  simulated.  These  shear-fracture  patterns  are  being  made  to 
help  understand  the  phenomenon  of  water  or  other  fluids  lubricating  an 
existing  fault  surface  under  large  lateral  loads.  The  calculations  made 
to  date  arc  discussed  later  in  this  report  and  apparently  do  not  agree 
with  previous  analytic  solutions  of  the  same  types  of  problems.  These 
differences  will  have  to  be  resolved.  Superposition  of  stress  fields 
along  with  diffusion  waves  in  the  strata  can  then  be  attempted. 


6 


RESULTS 


Several  calculations  have  been  made  with  the  code.  These  include 
simulations  o f  both  dynamic  ana  static  fractures.  The  stress  field 
around  a  static  cracK  (Fig.  1  and  2),  which  was  loaded  by  a  hydrostatic 
pressure  within  the  crack,  exhibits  concentrations  at  the  tip  of  the 
crack  and  at  both  sides  of  the  tip  concentrations  at  the  sides  of  the  tip 
are  areas  of  high  shear  stress.  The  medium  simulated  bad  tne  Lame’ 
constants  X  and  4  equal  to  10"  dynes  per  square  centimeter  and  Poisson’s 
ratio  equal  to  O.25.  This  calculation  was  performed  with  the  dynamic 
code  by  forcing  the  dynamics  to  settle  to  a  minimum  energy  conf iguration. 
Greater  detail  in  the  stress  field  can  be  expected  with  smaller  zones; 
however,  increasing  the  number  of  zones  will  also  increase  the  calcula¬ 
tion*  1  time.  The  c&lculational  mesh  used  had  15  x  19  grid  Intersections. 
Examination  of  the  contours  Indicates  that  the  stress  field  and  the 
distortion  are  not  symmetrical,  since  the  crack  was  not  placed  in  the 
center  of  the  medium,  A  similar  calculation  was  performed  with  the  crack 
discontinuity  located  symmetrically  and  a  symmetric  field  was  calculated. 

The  sawtooth  effect  seen  on  both  sides  of  the  crack  on  the  computer  plot 
is  caused  by  tne  routines  that  plot  the  contours,  not  by  the  dynamic  code. 
Other  contours  that  do  not  show  a  smooth  appearance  occur  in  areas  where 
the  contoured  function  has  a  shallow  gradient. 

Calculations  of  brittle  cracks  propagating  through  elastic  isotropic 
madia  have  been  performed.  A  simulation  of  a  tensile  crack  was  made  on  a 
grid  distorted  for  a  tensile  uniaxial  stress  field.  The  failure  criterion 
cnosen  for  crack  propagation  was  &  simple  tensile  break.  For  example,  the 
Iangrangian  zone  would  not  support  any  tensile  stress  perpendicular  to  the 
break  or  shear  stress  parallel  to  the  break  after  either  of  the  principal 
stresses  had  exceeded  a  certain  value.  Since  for  the  tensile  test,  the 
zone  would  not  show  a  crack  closure,  the  shear  stress  did  not  have  to  be 
modified  for  tangential  movement  of  the  crack  surface  in  the  event  that 
a  normal  force  of  compression  existed  within  the  zone.  The  uniaxial 
stress  in  the  medium  at  the  start  of  the  calculation  was  static  ,t  a 
value  of  1  kilobar.  Again  the  L*me;  constants  X  and  /a  were  10"  dynes/cm 
and  xoisson’s  ratio  was  0.25.  With  a  density  of  2.7  gm/cm  ,  the  congressional 


wav©  speed  was  0.3333  cm fp  sec,  and  the  distortional  wavs  ipeed  was 
0.1924  cm  I  p  sec.  The  start  of  the  crack  was  initiated  by  simulating 
the  situation  that  three  zones  along  one  boundary  of  the  grid  would 
suddenly  support  neither  tensile  nor  shear  stress  through  the  center 
of  the  zone  parallel  to  two  of  the  zone  sides.  Three  runs  were  made  with 
three  tensile  failure  levels  chosen,  1.3,  1.5  and  1.8  kilobars.  The 
crack  was  simulated  to  propagate  at  distortional  wav©  speed  for  the 
1.5  kllobar  tensile  strength  and  slightly  less  in  the  case  1.8  kilobar 
breaking  strength.  For  the  breaking  strength  of  1.3  kilobars,  the 
fracture  propagated  at  a  speed  faster  than  the  distortional  wave  speed. 

For  all  three  runs,  the  crack  moved  straight  across  the  grid  perpendicular 
to  the  uniaxial  stress  field,  showing  no  indication  of  branching  or 
splaying.  A  plot  of  the  elastic  flow  field  Is  shown  In  Figure  3  for  the 
tensile  breaking  strength  of  1.5  kilobars. 

Due  to  the  elastic  radiation  patterns  the  energy  released  from  tensile 
fracture  cannot  propagate  forward  of  the  crack  tip  at  speeds  greater  than 
the  distortion*!  wavs  speed.  The  crack  then  should  not  propagate  at 
speeds  greater  than  this,  either.  The  code  uses  a  finite  size  element, 
and  the  beam  affect  of  the  finite  calcul&tional  elements  may  be  a  factor 
in  the  hign  crack  velocity  for  the  smaller  failure  level  (1.3  kilobars). 

In  addition,  dispersion  could  exist  in  the  calculational  technique  because 
of  too  high  energy  release  rates  in  the  fracture  model,  which  may  also 
affect  the  fracture  velocity.  However,  the  failure  phenomenon  for  a 
dynamically  propagating  fracture  nay  not  be  so  simple  as  the  one  stated. 
Currently,  a  literature  search  is  under  way  to  help  resolve  this  and  other 
problems . 

As  a  preliminary  simulation  of  the  effect  of  water  lubricating  a 
sheared  medl/-ua  with  a  plane  of  weakness,  calculations  were  performed 
with  the  mesh  in  a  state  of  pure  snear.  The  medium  was  relaxed  with  an 
Imposed  shear  of  0.01,  corresponding  to  1  kllobar  shear  stress  using  the 
same  lause'  constants  and  Poisson's  ratio  as  before.  In  this  instance,.* 
plane  of  weakness  ms  chosen  laterally  through  the  mesh;  failure  could 
occur  only  along  thi  plane.  This  would  be  analogous  to  marking  a  pane 
of  glass  with  a  glass  cutter  and  Inducing  failure  along  the  scribed  line. 
The  failure  model  chosen  ms  such  that  when  the  shear  stress  along  this 
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line  exceeded  a  certain  value,  shear  fracture  would  occur  and  a  now 
"shear  stress  would  be  calculated  proportional  to  the  normal  force,  if 
compressive,  but  opposite,  to  the  iirection  of  motion.  In  addition, 
if  the  relative  motion  across  the  crack  approached  zero,  the  relative 
motion  was  stopped  and  a  static  frictional  force  was  used  for  the  shear 
force  in  the  zone . 

Simulations  on  the  sheared  medium  were  performed  with  this  failure. 

model.  The  fracture  was  started  by  removing  the  shear  force  in  three 

zones  on  one  boundary  of  the  mesh.  Three  calculations  were  made  with 

different  failure  levels  in  the  fracture  model.  The  levels  chosen  for 

the  parametric  study  were  1.3  T° _ ,  1.5  T  ,  and  1.8  T° _ _  where 

q  st  >>  1 s? 

T  is  the  initial  shear  stress  in  the  zone.  In  all  three  calculations, 

the  fracture  speed  stabilized  at  approximately  the  compressional  wave 

speed.  The  calculated  differences  in  the  three  simulations  caused  by 

the  changes  in  the  magnitude  of  fracture  or  slip  level  resulted  in 

changing  the  time,  required  to  build  up  to  this  speed. 

Examination  of  some  of  the  literature  indicates  that  the  shear 
fracture  should  not  propagate  at  the  compressional  wave  speed  and  that 
some  writers  [4]  propeje  a  fracture  toughness  term  as  the  fracture  speed 
increases.  A  term  of  the  form  AV^Txy/dx  was  includod  in  the 
fracture  model  to  simulate  the  resistance  to  fracture  as  the  fracture 
speed  increases.  When  this  term  was  included,  the  fracture  or  slip 
occurred  wiser  the  shear  stress  became  greater  than  a  minimum  fracture 
shear  and  greater  than  the  toughening  term.  In  addition,  fracture  always 
occurred  if  the  shear  stress  became  greater  than  some  maximum  level. 


In  equation  form,  fracture  would  occur  if 


and  Tn  > 

*y 


or 


> 


T 

max 


where  the  co-ordinate  X  is  taken  as  the  co-ordinate  of  the  craok  center 
lint  and  A  is  a  constant.  For  this  calculation  A  wa:  taken  as  0.08. 
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lining;  this  fracture  model,  the  fracture  speed  stabilized  at  about 
0.83  of  the  corapressional  wav©  speed.  This  calculation  showed  that  a 
more  complete  physical  description  of  dynamic  fracture  is  required. 

Contour  plots  of  the  shear  strain  indicated  that  it  did  not  Increase 
ahead  of  the  crack  tip  as  it  had  with  the  nonrate  dependent  term  included 
in  the  fracture  model  (Fig.  4).  Plots  of  the  velocity  vectors  of  the 
mass  points  provide  a  good  indication  of  the  elastic  radiation  patterns 
(Fig.  5).  These  patterns  for  shear  and  dilatational  radiation  are 
superimposed  on  the  plots,  but  close  examination  will  show  the  dilatation¬ 
al  patterns  quite  clearly.  A  magnification  of  the  grid  plot  (Fig.  6), 
which  is  actually  a  plot  of  the  medium  itself,  evidences  a  twisting  of 
the  crack,  along  with  a  tendency  to  open  slightly,  near  the  tip.  The 
stress  field  did  not  have  a  tensile  stress  imposed,  so  this  opening 
apparently  results  from  the  dynamic  crack  propagation. 

Other  calculations  included  simulations  of  slower  energy  release 
rates  from  the  crack  opening.  These  studies  were  performed  to  determine 
if  the  crack  was  releasing  energy  at  too  fast  a  rate  for  the  code  to 
calculate  properly.  Since  the  cod®  has  to  permit  a  crack  to  traverse 
an  entire  aone  at  the  instant  of  failure,  the  energy  release  rate  could 
be  causing  a  dispersion  of  the  energy  in  the  calculation.  These 
studies  of  slower  release  rates,  where  the  fracture  energy  was  released 
over  a  longer  period  of  time,  did  not  indicate  that  dispersion  distributed 
to  th*  calculated  high  crack  velocities.  The  effect  of  slower  energy 
release  rates  was,  again,  to  slow  the  rate  of  buildup  to  a  constant 
fracture  speed. 

All  the  crack  propagation  simulations  were  performed  with  a  30  x  60 
two-dimensional  mesh.  The  problems  were  set  up  so  that  the  fracture 
would  run  parallel  to  the  long  side  of  the  mesh.  With  the  fractures 
centered,  a  reflected  wave  from  a  parallel  side  of  the  mash  could  not 
impinge  on  the  fracture  discontinuity  until  84  microseconds  of  problem 
time  had  elapsed. 

Discussion  of  these  fracture  simulations  is  included  in  this  report 
for  completeness  and  to  illustrate  that  dynamic  fracture  is  a  complex 
phenomenon.  The  failure  criterion  used  here  apparently  does  not  describe 
the  phenomenon  of  dynamic  fracture. 
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B.  Cotterell  [4]  states  that  if  a  crack  is  propagating  in  an  infi¬ 
nite  plane,  the  crack  extension  force  will  continually  increase  unless 
the  crack  is  accelerating  all  the  tine.  Experimental  evidence,  although 
not  conclusive,  indicates  that  brittle  fractures  in  steel  quickly 
assume  a  steady  speed  of  propagation.  If  this  is  true,  the  fracture 
toughness  must  increase  with  the  stress  intensity.  It  is  unlikely  that 
such  a  phenomenon  could  continue  indefinitely;  at  some  stage,  the  excess 
sne.'gy  should  be  absorbed  by  a  fracture  branching. 

Preliminary  studies  of  the  literature  indicate  that,  analytically, 
the  study  of  moving  cracks  has  been  concentrated  on  two  basic  types: 

(a)  a  crack  of  constant  length  traversing  a  uniform  stress  field  at 
constant  velocity  [5]  »hd  (b)  a  crack  where  length  symmetrically 
increases  from  zero  with  constant  velocity  [4,  5]  •  Natural  cracks  are 
quite  obviously  of  type  (b),  since  typ®  (a)  assumes  that  the  fracture 
is  healing  beiiind  a  constant-length  crack. 

An  understanding  of  fracture  phenomena  in  geologic  strata  has 
fundamental  Importance  to  the  work  undertaken  in  this  contract.  The 
indications  of  the  analysis  perf oread  and  discussed  in  this  report  show 
that  a  complete  literature  search  is  Imperative  and  that  further  analysis 
and  correlations  must  be  made. 
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CCNCHfcICNS 


The  two-dimensional  elastic  dynamic  code  developed  for  use  in  the 
seismic  problems  produces  &  much  quieter  mesh  from  which  more  information 
can  be  extracted  than  do  some  of  the  existing  formulations.  Correlation 
with  the  source  detector  problem  (lamb's  problem)  on  an  isotropic  half 
space  shows  that  the  method  is  sound.  Simulations  of  fluid-driven 
fractures  and  a  fluid-lubricated  fault  surface  under  lateral  load  indicate 
that  single  fracture  models  do  not  provide  a  couplets  description  of  the 
dynamic  failure  phenomenon.  Further  development  work  on  fracture  models, 
with  an  accomp&ning  literature  search,  is  required. 
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Difference  Equations  For  Two-Dimeisional  Flow 
A.G.  Petschek  and  M.E.  Hanson 


ABSTRACT 

An  improved  technique  for  handling  two-dimensional  elastic  flow 
is  applied  to  several  physical  problems.  The  method  differs  from 
older  methods  in  that  no  unresisted  distortions  of  the  mesh  are  allowed. 
A  much  quieter  mesh  results,  and  more  information  can  be  extracted  from 
the  calculation. 
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APPENDIX 

FINITE  DIFFERENCE  EQUATIONS  FOR  THE  TWO  DIMENSIONAL 
ELASTIC  DYNAMICS  CODE  "TEMS" 

This  description  is  given  in  rectangular  coordinates. 

A.  THE  CALCULATIONAL  GRID  AND  MASS  ZONING 

A  two  dimensional  Langrangian  grid  is  placed  in  the  material^ 
dividing  the  material  into  rectangles.  Figure  7  is  a  schematic  of  the 
grid  and  a  numbering  system  for  the  center  and  corners  of  the  zones 
used  in  this  description  of  the  code.  In  particular 
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The  mass  at  each  vertex  of  the  rectangles  is  calculated  initially  (at 
time  zero)  and  held  constant  for  the  entire  calculation.  This  ensures 
conservation  of  mass  in  the  calculation.  The  mass  concentrated  at 
each  vertex  is  taken  as  1/4  the  sum  of  the  masses  of  the  adjoining 
zones.  For  example: 

m©=  '( Da® 

m,  ,  =  —  (m_  +  m _ +  m_  +  m  ) 

w  4  ©  ©  ®  © 

The  masses  at  the  other  vertices  are  calculate*  similarly. 

The  area  of  a  quadrangle  is  taken  as 


(Al) 
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A&  %* 1 V® 


...n  1  r  n.  n  n.  n,  n  n.  n.n  n.l 

(AI 2  tx2  y3  ‘  y4  1  +  *3  <V4  '  y2  >  +  *4  (y2  ‘  yJ  >J 


(A2) 


.n  1  r  n.  n  n%,  n.  n  n.  ,  n,  n  n.l 
(AuVt3=  2  lX2  (yi  ■  yl  >  +X4  (yl  -  y2  >  +X1  y2  '  y4  >J 


Aj  and  are  the  areas  of  the  triangles  I  and  II. 


B.  STATE  EQUATIONS 


The  compression  at  cycle  n  is  given  by 


n  A  0 

n  _  jp _  _  A_ 

p  0  i  n  i 
r  .A 


(A3) 


The  code  permits  the  effective  bulk  modulus  to  be  given  by  an  expansion 
of  the  form 

k"  =  =  a  +  2b  (l).n  -  1)  +  3c  (1)"  -  l)2  +  4d  (1)."  -  l)3  (A4) 

where  the  coefficients  b,  c,  d  . . .  are  empirical  fits  to  experimental 
shock  data  and  can  be  found  in  the  literature,  for  example  Walsh  et  al. 

[7]  .  The  coefficient  a  is  the  linear  bulk  modulus  of  the  material. 

Given  one  of  the  Lame  constants  and  Poisson's  ratio  the  remaining  Lame 
constants  can  be  found  for  linear  elasticity.  Poisson's  ratio  is  taken 
to  be  constant  so  that  the  Lam£  "constants"  are  given  by 


n_  3  T,  n  ,  1  -  2  v 

=  2  1  {mr 


(A5) 


x."  .  2£iJ 
1  -  2 


where  fi  is  the  shear  modulus  and  v  is  Poisson's  ratio. 
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C.  STRAINS  AND  STRESSES 


Consider  a  rectangle  whose  original  corners  were  at  (X,  Y), 

(X  +  h,  Y),  (X  +  L,  Y  +  K),  (X,  Y  +  K)  and  let  the  corners  presently 
be  at  (x,  y),  (x  +  L^,  y  +  k^),  (x  +  Ly  y  +k^).  (x  +  L^,  y  +  k^)  as 
6hown  on  Figure  8  . 

Now  consider  a  point  whose  original  coordinate  was  (X  +  U, 

Y  +  V)  and  whose  present  coordinate  is  (x  +  u,  y  +v).  Expand  u  and 
v  as  in  Eq.  (A25). 

The  coefficients  a  through  f  can  be  obtained  by  solving  Eq.  (A26) 
which  gives 

&n  =  l-  /L  bn  =  lA  /K 

C  4 

dn  =  k2/L  en  =  k4/K  (A6) 

°n  =  ( *3  “  *2  "  V/LK  =  {k3  '  k2  "  V/LK 

The  distortion  of  the  grid  will  produce,  in  general,  both  a 
rotation  and  a  distortion  of  each  individual  mesh.  As  the  rotation 
will  lead  to  strains  which,  by  Hooke's  law,Eq.  (A10), produce  large 
stresses,  it  is  necessary  to  account  for  the  rotations  separately.  In 
the  distorted  mesh,  a  line  of  constant  V  has  a  slope  given  by  Eq.  (A27). 
Rotating  the  distorted  mesh  through  the  angle  -Q  putt  this  line  parallel 
to  its  original  orientation.  After  rotation  through  this  angle,  the  ooint 
originally  at  U,  V,  is  at  u',  v'  where 

u'  =  u  cos  B  +  v  sin  B  (AT) 

v'  =  v  cos  B  -  u  sin  B 
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The  strains  are  accordingly 


U  ?)n  =  ~a(-~irr^~  =  (a  +  cV)ncos  0n  +  (d  4  fvfsin  fl"  -  1 

x'  |J  u 


Uy?)n  =  •a(yl^-vV)-  =  (e  +  fU)ncos  -  (b  +  cU)nsin  0n  -  1 


(A8) 


(ry,f.)n  -  a(V--^V)n  +  -  (b  +  cU)ncos  0  "  +  <e  +  fU>\in*  “ 


Eliminating  the  trigonometric  functions  by  use  of  Eq.  (AZ7)  gives 


(  \/(a  +  cV)Z  +  (d  +  fV)Z  -  1)" 

~(ae  -  bd)  4-  U(af  -  dc)  4-  V(ec  -  bf)  _  A  1 
\/(a  +  cV)2  +  (d  +  fV)Z 

(b  t  cU)  (a  4-° cV)  4-  (d  4-  fV)  (e  4-  fU)*[  " 

/  2  Z  I 

L  V(a+cV)  4  {d  4-  fV)  J 


(A9) 


using  Hooke's  Law,  the  stresses  at  a  point  along  the  line  of  constant 

V  are 


,v/  =  ^+MB(«y?)UVx?  >» 

.  .n  n.  Q  .n 
(<y'x')  ^y'x' 


(A10) 


Similarly  a  line  of  constant  U  in  the  distorted  mesh  has  a  slope 
given  by  Eq.(A28).  Rotating  the  distorted  mesh  through  the  angle 
-^puts  this  line  parallel  to  its  original  orientation.  After  rotation 
the  point  origins  .ly  at  U,  V  is  at  u',  v'  where 
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(All) 


u'  =  u  cos  <p  -  v  sin 
v'  -  v  cos  <f>  +  u  sin 


the  strains  arc  accordingly 


(€  f  )n  =  ■a(uyuU)-  =  (a  +  cV)n  cos  <£n  -  (d  +  fV)n  sin  $  n  -  1 

(€  ^  )n  =  —  ^  ~  •  =  (c  +  fU)ncos  +  (b  +  cU)nsin<£n  -  1  (A12) 

y '  d  * 

(rxi.  ,"  =  iLv'au^+M"Md  +  fv)ncosr 


+  (a  +  cV)nsin^“ 


Again  eliminating  the  trigonometric  functions  with  Eq.  (A28)  gives 


.  6  ,n  T(ae  -  bd)  +  U(af  -  dc)  +  V(ec  -  bf)  , 

ux.  )•  ■  —  =  Z _  _  "  1 

*/{e  +  fU)2  +  (b  +  cU)2 


<€  t  )"  =  (V£7fu)2  +  (b  +  cU)2  -  l)n 


i 

A  %n_  f(a  +  cV)  (b  +  cU)  +  (d  +  fV)  (c  +1 

<yxy  -  - - - 

V(e  +  fU)  +  (b  +  cU) 


(A13) 


The  stresses  at  a  point  along  the  line  U  =  constant  are 


K  =  t2/i  +  X)n(«  t  )"  +  )n 


<”xv>n  *  f{rxt  >" 


(A14) 


The  angle  is  taken  positive  in  the  clockwise  direction  and  the  angle 
0  is  taken  positive  in  the  counterclockwise  direction. 
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Occasionally  it  has  been  desirable  to  recover  the  stress  field 
in  the  original  frame  in  which  case  the  stresses  in  the  interior  of  the 
zone  arc  defined  by 

(°xx*n  =  (<x'x',n  COS  ^  +Kf-.y.)n  8in2  ^  "  2  (Txiy»)n  8in^  co®  t 

(O'yy)”  =  (O’y.y.)”  C0S  ^  +  i^x,)"  ^  ^  +  2  (  ^  ,  )"  8™  »J  COS  ^ 


(A15) 


where 


Tx|yi  E  (c  ^  ^x'y'^  ^ 


tan^  =  (tan0n  -  tan  ^  'j/Z 


(A16) 


This  amounts  to  averaging  the  rotation  and  shear  stresses 
computed  along  the  lines  U  =  constant  and  V  =  constant. 

D.  THE  MOMENTUM  EQUATIONS 

For  the  dynamics,  the  stresses  are  integrated  on  the  quarter 
zones  about  the  mass  point  k,  L  .  For  example,  the  integration  is 
carried  out  along  the  dashed  line  segments  about  the  mass  point  k,  L 
in  Figure  7. 

We  must  integrate  the  appropriate  stresses  along  the  line 
V  =  K/2  over  the  segment  0<U<L/2  or  L/2<U<L  depending  on 
the  zone  corner  bordering  on  the  mass  point.  Likewise  the  integration 
must  be  carried  out  for  the  appropriate  stresses  along  the  Jfir.o  U  =  L/2 
over  the  segment  0<V<K/2  or  K/2<V<K.  Either  the  strain  defi  .itions 
(A8)  and  (A12)  or  (A9)  and(Al  3)  maybe  used.  We  will  use  the  definitions 
(A8)  and  (A12)  to  correspond  to  the  form  in  the  text  of  the  report.  The 
integrations  in  cither  case  arc  trivial 
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(A17) 


(Fy'y'lin=/2  [<2M  +X){<e  +  fU)  cos  8  -  (b,  +  cU)  sin  Q  -  }  V 
X1 


_  •  • 

+  X|(a  +  cK/2)  co  8  +  (d  +  fK/2)  sin0  -  lj  ]  dU 


(Fy^,)"  =  /  (b  +  cU)  cos  8  +  (e  +  fU)  sin  0}]  dU 
X1 


where  x^  and  are  the  appropriate  limits  either  x^  =  0,  x^  =  L/2 
or  Xj  =  L/2,  x^  =  L,  and 


X4 

(Fx,xI)in  =  /  [(2/i+  X)  /(a  +  cV)  cos  (f>  -  (d  +  £V)  sin  -  l} 

X3 


n 

+  X|(c  +  fL/2)  cos  +  (b  +  cL/2)  sin  ^  -  l|J  dV 


(Alb) 


x  n 

Rd  +  £V)  cos  $  +  (a  +  cV)  sin  d>|J  dV 


where  x^  =  0,  =  K/2  or  =  K/2,  =  K  depending  on  the  position 

of  the  ?,one  about  the  mass  point.  Then 

n 

(Fy|y^.n  =  (2/t  +  X).n[(e  +  afL/4)  cos  8  -  (b  +  dcL/4)  sin0  -  lj  L/2 

n 

+  X.n  [(a  +  cK/Z)  cos  8  +  (d  +  fK/2)  sin  8  -  l]  L/2  (A19) 


n 

(F  ^,).n  =  fi  ,n  |(b  +  OcL/4)  coe  0  +  (c  +  afL/4)  sinOj  L/2 
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and 


(F  »  (2/l+X).n  [(a  tacK/4)  cos  <f>  -  (d  +a  fK/4)  sin  $  - 

X  X  l  l  L 


n 

ll  K/Z 


n 

+  X-  £(e  +  fJj/Z)  cos  &  +  (b  +  cL/2)  sin  <j>  -  lj  K/2  (A20) 


n 

(F  ,  ,).n  =  /i  [(d +afK/4)  cos<£  +(a  +acK/4)  sin<£  ]  K/2 

x  y  1  i  l  J 


where  the  subscripts  i  denote  the  zone  centers,  for  example  those 
numbers  enclosed  in  circles  in  Figure  7,  in  relation  to  mass  point 
k ,  L  ,  and  n  is  the  cycle  number. 


a 

a 

a 

a 


=  1  for  Equations  (A19)  for  zones 
=  1  for  Equations  (A20)  for  zones 
=  3  for  Equations  ( A 1 9 )  for  zones 
=  3  for  Equations  (A20)  for  zones 


and 

and 

and 


As  in  the.  text,  the  first  subscript  designates  the  normal  to  the  surface 
across  which  a  stress  in  the  direction  of  the  second  subscript  is 
exerted.  The  primes  on  the  subscripts  denote  that  we  are  in  a  frame 
rotated  with  respect  to  the  original  frame. 

We  must  then  rotate  these  forces  back  to  the  original  frame  to 
apply  them  to  the  momentum  equation. 


(F  ,  )n  =  (F  ,  ,)n  cos  9  n  -  (F  ,  ,)n  sin  9  n 
y'x  l  y'x'  l  y'y*  l 


(F  ).  =  (F  ,  ).  cos  9  +  (F  ,  ,).n  sin  9 

y  y 1  y  y*  y'x'  1 


(A21) 
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and 


(F  ,  )"=(F  ,  ,)"  cos  V>"  +  (F-  ,  ).n  sin<f>  n 

'  x'x  x  x'x'  i  r  x'y'  1  T 


(F  ,  )."  =  (F  ,  cos  <f>n  -  (F  ,  ,).n  sin  p  ” 
x'y  1  x'y'  1  r  x'x'  x  r 


(A22) 


The  velocities  can  be  computed  from  the  force  field  by 


n  +  1/2  .  n  -  1/2 

X  r  Y 

k  A  k,A 


+  f(F*.*V-  <F*.xK‘  (F,XU 
1  xx©  xx©  x  ©  © 


n 

+  (F  .  W(F  ,  )  -  (F  .  )  -  (F  .  )  "I  Atn/m,  1 

y'xQ  y'x0  y  x0  y  X@J  k*' 


.  n  +  1/2  n  -  1/2 


*  Va"  '  *  4  K'y>®+  (Fy'y&-  (Fy'y&  (Fy'y>© 
+  (Fx'v^T)’  (Fx'y^)’  (Fx'y^j)+ 


The  positions  can  be  incremented  by 


n  +  1  n  ,  n  +1/2.  n  +  1  /2 

X'k.A  Xk,  L  4Xk,.  At 


n  +  1  n  n  +  1/2.  n  +  1/2 

yk,t  =  Fk.A  +  Vi  A‘ 


(A23) 


(A24) 


E.  TIME  STEP 

In  accord  with  the  Courant  criterion,  the  time  step  is  taken  to 
be  the  m*nimum  over  the  mesh  of 
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At 


n  +  1  /2 


=  CAr 


n  +  1 


/a 


where 


Ar 


n  +  1 


An  +  1 


S  is  the  longest  diagonal  of  the  quadrangle,  a  is  the  local 
dilitational  wave  speed  and  C  has  been  taken  to  be  0.  6  . 

F.  THE  MAPPING  TECHNIQUE 


Evidently,  one  ought  to  difference  the  equations  so  that  there 
are  no  unresisted  distortions  of  the  mesh.  This  has  been  done  as 
follows:  consider  a  quadrilateral  which  was  originally  rectangular 
and  had  its  lower  left  hand  corner  at  X,  Y.  Let  the  lower  left  hand 
corner  now  be  displaced  to  x,  y  and  the  point  originally  at  X  +  U, 

Y  +  V  be  displaced  to  x  +  u,  y  +  v.  Expand  the  present  position  as 
a  function  of  the  original  position  as  follows 


u  =  aU  +  bV  +  cUV  +  . . . 
v  +  dU  +  eV  +fUV  +  .  .  . 


(A25) 


Just  six  parameters  in  the  expansion  have  been  written  down  because 

this  is  enough  to  specify  completely  the  positions  of  the  corners  and  is 

as  many  as  can  be  determined  by  those  positions.  The  term  proportional 

2  2 

to  UV  is  preferable  to  one  proportional  to  U  or  V  because  it  maps 
the  edges  of  the  original  rectangle  into  straight  lines  and  thir  ensures 
that  the  meshes  of  the  distorted  lattice  fit  together  without  gaps  or 
overlap. 

The  coefficients  a  through  f  may  be  found  from  the  coordinates 
of  the  grid.  That  is,  if  the  intersections  of  grid  lines  around  the  mesh 
are  numbered  from  1  to  4  as  before,  the  one  with  U  =  V  =  0  is  numbered 
1,  and  the  original  rectangle  had  the  dimensions  L  by  K,  then  a 
through  f  are  found  from 
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u2  =  aL 


v2  *  dl 

u  =  aL  +  bK  -f  cKL 

3  (A26) 

-  dL  +  eK  +  fKL 

u.  =  bK 

4 

v .  =  eK 

4 

Given  a  through  f  ,  one  may  compute  the  dieplacments  u  -  U 
and  v  -  V  and  from  tne  displacements  the  strains  as  a  function  of  position 
throughout  the  mesh. 

To  compute  accelerations,  the  mesh  is  bisected  in  the  X  and  Y 
directions,  and  the  masses  of  the  four  adjacent  quarter  zones  as  associated 
with  each  mesh  point.  The  force  is  comouted  by  integrating  the  stress 
around  the  circumference  of  the  mass.  To  do  this  the  stresses  on  the  line 
segments  U  =  L/2,  0<V<K/2;  U  =  L/2,  K/2<V<K  and  so  forth  are 
required.  These  can  easily  be  computed  from  the  strains,  provided 
proper  account  is  taken  of  the  rotation  of  the  mesh.  To  do  the  latter, 
we  rotated  the  zone  until  the  line  along  which  the  stress  was  to  be  com¬ 
puted,  i.  e.  V  =  K/2  or  U  =  L/2  was  parallel  to  its  original  direction, 
i.  e.  the  X  or  Y  axis.  This  implies  a  rotation  through  an  angle 


$  =  tan 


-1  d  +  fK/2 
a  +cK?2 


(A27) 


in  the  former  case  and  an  angle 


9 


=  tan 


-1  b  +  cL/2 
e  +  fL/2 


(A28) 


in  the  latter.  If  the  rotations  are  not  taken  into  account,  Hooke's  law 
as  used  below  will  give  rise  to  large  fictitious  stresses. 
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After  rotation  the  strain  in  the  x  direction  along  the  line  of 
constant  V  is 


9 


(a  +  7  cK)  cos  9  +  (d  +  7-  fK)  sin  9-1. 


<A29) 


The  primes  indicate  that  a  rotated  coordinate  system  is  in  use. 

The  stresses  are,  in  the  case  of  plane  strain  (no  distortion 
perpendicular  to  the  plane  considered)  by  Hooke's  law,  for  instance 


cr  ,  ,  =  (2  /i  +  X  )  £  (e  +  fU)  cos  9  -  (b  +  cU)  sin  9  -  lj  -*■ 
X  [(a  +  j  cK)  cos#  +  (d  +jfK)  sind  -  l] 


(A30) 


Where /i  and  X  are  Lame  constants  and  the  first  subscript  designates 
the  normal  to  the  surface  across  which  a  stress  in  the  direction  of  the 
second  subscript  is  exerted. 

The  corresponding  forces  are  obtained  by  integrating  along  the 
appropriate  half  bisector  of  the  mesh  to  give,  for  instance 


y'y' 


J  dlT  <r  ,  ,  =  (2/a+  X )  £  (e  +  ^  faL)  cos  9  -  (b  foL)  sin  9  -  lj 
+  X[(a  +  jcK\  cos  9  +  (d  +  |-fK)  sin  9  -  1  ]  ~L 


(A31) 


where  the  limits  on  the  integral  are  0  and  L/2  or  L/2  and  L,  and 
a  is  1  in  the  former  case  and  3  in  the  latter.  Then  the  force  is  rotated 
back  into  the  original  Eulenan  frame,  that  is 


=  F  ,  , 

cos  9  +  F  . 

sin  9 

y'y 

y'y' 

y'x1 

=  F  ,  , 

sin#  +  F  . 

.  cos  9 ' 

y'x 

y'y' 

y'x' 

(A32) 
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Finally  the  force  is  summed  over  the  eight  line  segments  surrounding 
the  intersection,  and  is  used  to  accelerate  the  mass  centered  there. 

G.  GENERAL 

The  boundaries  of  the  mesh  can  be  treated  as  having  phantom 
zones  with  no  tractions,  normal  forces  or  masses.  The  code  will  not 
correctly  calculate  displacement  disturbances  with  frequencies  shorter 
than  the  time  acquired  for  a  dilatational  disturbance  to  cross  4  zones. 

If  these  disturbances  are  produced  on  the  boundaries  of  the  mesh,  dis¬ 
persion  of  the  induced  wave  trains  is  observed. 

A  tensor  viscosity  has  been  developed  and  used  with  the  code. 
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Firure  1  Grid  for  static  crack  with  distortion  magnified  by  10.  Hydrostatic  pressor® 
loo  the  crack  was  taken  as  1  kilobar,  lame'  constants  for  media  were 
X  =  o  =  0.1  meganars. 


•.00 


0 


*  00 


s  c  i  g  k  s  8  e 


tiOG  budd 


!  R  S  rt  IB  K  «  K  C. « 

aiGOUJlW  966SV0nfiI  3WI1 


9  £22  313JL3 


Figure  2  Contours  of  P-Q  for  the  internally  loaded  static  crack 
the  principal  stresses.  Contour  interval  is  125  bars, 
125  bars,  maximum  contour  is  3500  bars,  the  'sawtooth  a 
of  the  crack  is  due  to  the  plotting  routines.  No  init 
on  the  medium.  Contours  are  of  the  *rid  snown  on  fiau 
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Figure  3  Vector  velocity  plot  of  mass  points  or  elastic  flow  field  around  tensile 
crack  in  a  uniaxial  strained  media.  Time  slice  of  plot  is  28.8  p  sec 
after  crack  was  initiated.  Vector  len-’-th  magnification  is  625.  ^one 
size  for  calculation  was  1  centimeter  square.  Lame*  constants  for  the 
media  twre  A  =  p  =  0. 1  megabars. 
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Figure  5  Velocity  vector  plot  or  elastic  flow  field  around  the  shear  fracture  in 
the  media  under  pure  shear  when  crack  wp.s  initiated.  Time  slice  war 
9^.^  v  seconds  after  crack  was  initiated.  Fracture  model  contained 
toughening  term  described  in  the  report.  Vector  length  magnif ication  is 
250.  Inme'  constants  for  the  media  were  X  =  y  =  0.1  megabars 


Figure  6  Grin  plot  of  media  during  shear  fracture  in  the  media  under  pare  shear. 

Time  slice  was  u  seconds  after  crack  was  initiated.  All  distortions 

were  magnified  by  5^-  Note  twisting  effect  at  crack  tip. 
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Figure  8.  Comparison  of  undistorted  and  distorted  meshes  indicating  the 
meaning  of  u  and  v  . 


